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We numerically solve the Boltzmann equation for trapped fermions in the normal phase using 
the test-particle method. After discussing a couple of tests in order to estimate the reliability of 
the method, we apply it to the description of collective modes in a spherical harmonic trap. The 
numerical results are compared with those obtained previously by taking moments of the Boltzmann 
equation. We find that the general shape of the response function is very similar in both methods, 
but the relaxation time obtained from the simulation is significantly longer than that predicted by 
the method of moments. It is shown that the result of the method of moments can be corrected by 
including fourth-order moments in addition to the usual second-order ones and that this method 
agrees very well with our numerical simulations. 
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I. INTRODUCTION 

In experiments on ultracold trapped Fermi gases, there 
are many situations where the system is out of thermal 
equilibrium. The first one is of course the trapping and 
cooling stage, i.e., before the system has reached its equi- 
librium state which is usually the starting point for the 
actual experiment. Then, in some experiments the sys- 
tem is excited in order to observe its dynamical behavior. 
For instance, many experiments studied collective oscilla- 
tions of the system PHD] , another example being a recent 
experiment at MIT where the collision of two atom clouds 
(both in equilibrium) was studied Finally, often the 
system is not imaged directly during the experiment, but 
only after the trap was switched off and the system has 
expanded for a certain time, in order to increase its size. 

The modeling of such time-dependent processes from 
the theoretical point of view can be quite complicated. 
For practical reasons, only semiclassical approaches are 
suitable for the description of time-dependent phenom- 
ena involving typically several 10 5 atoms in a three- 
dimensional, non-uniform geometry. In some cases, it 
is possible to use hydrodynamic approaches: Supcrfluid 
hydrodynamics describes the expansion 0] and the col- 
lective modes of superfiuid systems at zero temper- 
ature. Hydrodynamics is also applicable in the normal- 
fluid phase if the mean time between collisions is much 
shorter than all other time scales of the process under 
consideration, so that the system can always be consid- 
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ered to be in a local equilibrium (lfj, . Superfiuid and 
normal hydrodynamics can be combined to two-fluid hy- 
drodynamics in order to describe superfiuid systems at 
finite temperature [l2j . However, in many cases hydro- 
dynamic approaches are not sufficient. In all cases where 
it is important that even locally the distribution /(r, p, t) 
of the atoms is not an equilibrium one, the Boltzmann 
equation allows a very general description, provided the 
system is in the normal phase. If the system is supcr- 
fluid, a more elaborate theory is necessary which couples 
the dynamics of the quasiparticle distribution function to 
the dynamics of the superfiuid order parameter [ill [l4| ■ 

In the past, several authors used the Boltzmann equa- 
tion for the investigation of collective oscillations in 
normal-fluid trapped Fermi gases H EH EMU- In most 
cases, the Boltzmann equation was not solved directly 
in order to find the distribution function /, but semi- 
analytical approximate solutions were found by u sing the 
scaling ansatz 11 1 or the method of moments [fj|. ll7H19| . 
These methods rely explicitly or implicitly on the as- 
sumption that the collision term can be treated in the 
relaxation-time approximation, with a single relaxation 
time t which is independent of the position in the trap. 
An exception is the work by Toschi et al. [HI, EH , where 
the Boltzmann equation was solved numerically, using a 
test-particle method very similar to the one we are us- 
ing here. The test-particle method for the solution of 
the Boltzmann equation has been used for many years in 
nuclear physics for the simulation of heavy-ion collisions 
[20j . In the context of trapped atoms, it was also used 
for the simulation of the dynamics of the thermal cloud 
in a Bose-Einstein condensate [2l[ and (without colli- 
sion term) of the normal-fluid component in a supcrfluid 
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Fermi gas fl4| . 

For the collision term of the Boltzmann equation, it 
is important to know the cross section which in princi- 
ple can be modficd by in-medium effects. In the work 
by Ricdl ct al. it was shown that by using in the 
Boltzmann equation instead of the free cross-section the 
in-medium one, the agreement between theoretical and 
experimental frequencies and damping rates of different 
collective modes is deteriorated. The reason is that the 
in-medium cross section is larger than the free one, so 
that the relaxation time is dramatically reduced. In our 
previous work [Is}, our aim was to include in addition 
to the in-medium cross section medium effects into the 
mean-field potential. This mean field resulted in better 
density profiles and allowed us to understand the shift 
of the quadrupole frequency in the collisionless regime at 
very low temperature observed in [J]. However, it did not 
help to improve the agreement between theory and ex- 
periment in the region of higher temperatures, where the 
properties of the collective modes are completely domi- 
nated by collisions. 

Hence, one of our motivations for the present work was 
to check the validity of the relaxation-time approxima- 
tion which is implicitly made in the method of moments. 
Here we will restrict ourselves to the case without mean 
field and with the free cross section. As we will show in 
Sec. IIII1 the numerical solution of the Boltzmann equa- 
tion gives indeed a significantly longer relaxation time 
than the method of moments. As we will show, this 
discrepancy is due to the restriction of the method of 
moments to second-order moments in the existing litera- 
ture [tl Il7l - fl9j . Once fourth-order moments are included, 
the results of the method of moments and of the numer- 
ical solution are in good agreement. However, already 
in the simplest case of a spherical harmonic trap with- 
out mean field, the inclusion of fourth-order moments is 
a very tedious task, while the numerical method can be 
generalized to more realistic cases. 

In addition, there are some other reasons why we felt 
the necessity for a numerical method. For example, there 
are damping effects due to the anharmonicity of the trap 
potential which cannot be described by the method of 
moments. Another advantage of the numerical method 
is that it offers the possibility to simulate not only the os- 
cillation of the cloud, but also the subsequent expansion 
after the trap has been switched off. 

In Sec. |TT] of the present paper, we give a detailed de- 
scription of the method. In particular, we explain in 
detail how the collisions are simulated, since our method 
is somewhat different from that of Ref. [HI]. Moreover 
we discuss some tests we made in order to estimate up 
to which precision we can trust our simulation. Then, 
in Sec. IIII1 we come to the main point of our article and 
calculate the properties of some collective modes for a 
system in a spherical harmonic trap. While the sloshing 
and breathing modes are rather trivial, the frequency and 
damping rate of the quadrupole mode are very sensitive 
to the collisions. We compare the numerical results with 



those of the method of moments. Finally, in Sec. IIVI we 
summarize and give an outlook to future studies. 

Throughout the paper, we will use units with h = 
ttB = 1 (h = reduced Planck constant, ks = Boltz- 
mann constant). The strength of the interaction is char- 
acterized by the dimensionlcss quantity kpa 1 where a is 
the scattering length. Concerning the Fermi momentum 
kp and the Fermi energy Ep we follow the usual con- 
vention that these quantities are defined by the corre- 
sponding ones of an ideal Fermi gas at zero tempera- 
ture, i.e., kp = \/2mEp, m being the atomic mass, and 
Ep = (3A) 1 / 3 wo, where N is the number of atoms and 
wo the trap frequency. Temperatures will be measured in 
units of the Fermi temperature Tp = Ep (since ks = 1). 



II. DESCRIPTION OF THE NUMERICAL 
METHOD 

A. Test-particle method 

We study a two-component (a =f, I) gas of fermionic 
atoms of mass m in a potential V(r, t) with attractive 
interaction a < 0. We assume that the system is in 
the normal phase and that it can be described semiclas- 
sically by phase-space distribution functions f a (r,p,t). 
In this paper, we will restrict ourselves to the case that 
the distribution functions of both spin states are equal 
(/f = fi = /), but the generalization of the method to 
the cases of different distribution functions, more than 
two components, or components with different masses is 
straight-forward. 

The time evolution of the distribution function / is 
governed by the Boltzmann equation [22[ 

./ + !•• V r / + p-V p / = -/[/], (1) 

where the left-hand side (lhs) describes the particle prop- 
agation, with 

r=— and p = -W, (2) 

TO 

and /[/] on the right-hand side (rhs) denotes the collision 
term which will be discussed later. The potential felt by 
the particles is the trap potential that contains a static 
part and a time dependent one (which will be used to 
simulate the excitation of the collective modes) V(r, t) = 
V T (r) + Vi(r,t). 

The density per spin state is related to the distribution 
function by 

p(r,t)= J -^f(r,p,t), (3) 
and the number of atoms is given by 

N = N t + N l = 2 [ d 3 rp(r,t). (4) 
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The basic idea of the test-particle method (also called 
pseudopartiele method) for solving the Boltzmann equa- 
tion consists in replacing the continuous distribution 
function by a sum of delta functions, 



/(r, p,t) = ^ ]T(2^) 3 <Kp - P 4 (*))<*(r - *(«)) , (5) 

i— 1 

where N is the number of "test particles" . This allows 
one to express the average of an arbitrary single-particle 
observable F(r,p) in the simple form 

2 f d 3 rd 3 p , . „, . 1 . 

{f) = nJ ^^^ f ^-jT, f ^)- 

(6) 

In order to sample the six-dimensional phase space, it 
is necessary to choose a sufficiently large number of test 
particles N (usually N > TV). Neglecting the collision 
term /[/] for the moment, it is easy to see that Eq. §5§ 
satisfies the Boltzmann equation (Eq. (fT])) if the posi- 
tions Yi and momenta p i of each test particle i follow the 
classical equations of motion, Eq. ([2]). 

In practice, the delta functions in Eq. ([5]) can pose 
some problems. For instance, they do not result in a 
continuous density p(r). Therefore it is often useful to 
replace them by Gaussians of width w r and w p in position 
and momentum space, respectively: 



N 



S(jP - P»)tf(r - n) -> g w (p - Pi)g Wr (r - n) , (7) 



with 



The widths w r and w p must be adapted such that they 
smooth out the fluctuations due to the finite number of 
test particles, but not the structure of the distribution 
function /. The statistical fluctuations are of the order 
of ^Nwf.Wp/N) -1 / 2 , i.e., the first condition is equivalent 
to 



w r w p 3> 



N 
2N 



1/3 



(9) 



The second condition implies of course that w r <C Rtf 
and w p «pf, where Rtf and pf are the Thomas-Fermi 
radius and the Fermi momentum, respectively, but this 
is not always sufficient. At low-temperature, it is crucial 
to resolve the rapid change of the distribution function 
around the Fermi surface, i.e., 



w p < pf 



T 
T~f 



and w r <C Rtf 



T 
T~f 



(10) 



In practice, as the computation time increases as N 2 , it 
turns out that the conditions © and (fTOj) cannot simul- 
taneously be satisfied at too low temperatures. 



B. Particle propagation 

In the absence of collisions, the numerical task consists 
only in solving simultaneously the classical equations of 
motion ^ for the N test particles. We do this by using 
the velocity Verlet algorithm J24|, which contrary to the 
original Verlet algorithm [25[ uses the positions Yi(t n ) 
and velocities ~Vi{t n ) = p^t^/m as starting point for the 
time step from t n to t„+i — t n + At. The propagation 
from t n to t n+ i is done according to 

Vi(i„+i/ 2 ) = Vi(tn) + aa{t n )At/2 (11) 
ri{t n+1 ) = n(tn) + v l (t n+1/2 )At (12) 
Vi(* n +i) = v i (t n+1/2 ) + ai(i n+ i)At/2 , (13) 

where aj(t) = —'W(ri(t),i)/m is the acceleration of the 
i-th test particle. If it is written in this way, it is ob- 
vious that the velocity Verlet algorithm is identical to 
the leap-frog algorithm [26j]. Note that the accelerations 
aj(t„ + i) can be reused in the next time step, so that the 
algorithm needs only one evaluation of the acceleration 
per time step, but nevertheless its global error is of the 
order O(At) 2 . This allows us to obtain a good accuracy 
for reasonable time steps At. A good test of the particle 
propagation is to check the energy conservation: typi- 
cally we find \E t (t) - Ei(0)\/Ei(0) ~ 1CT 6 for each test 
particle and for all times considered. 



C. Collision term 

The rhs of the Boltzmann equation (TTJ describes the 
collisions between particles of opposite spin. It thus de- 
pends on the scattering cross section da/dQ and reads 

/|/l ^/ll/ <,f! S |v - v ' l[// ' ,1 - / ' )(1 - /:) 
-rawxwi)]. (M) 

In the first term, p and p 1 are the incoming momenta, p' 
and p^ are the outgoing ones, f2 is the solid angle formed 
by the incoming relative momentum p — p 1 and the out- 
going relative momentum p' — p[, and / = /(r, p, t), 
fx = /(r, Pi,i)i etc. In the second term, the role of in- 
coming and outgoing momenta is exchanged. Momentum 
and energy conservation implies p + p 1 = p' + p[ and 
|p — p-J = |p' — p[\. Here we consider the case of pure 
s-wave scattering, in which the cross section is isotropic, 
i.e., da/dQ = a/Air. In principle the cross section is mod- 
ified by medium effects [3 HH , but in the present paper 
we will only use the free cross-section (i.e., the cross- 
section for the scattering of two atoms of opposite spin 
in free space) which is given by (23j 
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1 + (qa) 2 ' 



(15) 
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where g = |p - Pl |/2 =|p'-pi|/2. 

In our numerical simulation, the collision term is in- 
cluded by allowing the test particles to collide with each 
other. The cross section of the test particles, a, is re- 
lated to the cross section of the atoms by a = aN/2N 
(since N test particles represent N/2 atoms of a given 
spin). Whether a pair i,j of test particles collides in 
a time step t n or not is determined as follows: First, 
we determine if the two particles are at their closest ap- 
proach in the present time step. Explicitly, if we write 
r,jj = r, — Yj and Vy = Vj — Vj, the closest approach 
is reached at t m i„ = t n — • vy/ify- and we check if 
\tmin — t n \ < At/2. If yes, we calculate the correspond- 
ing minimal distance by d 2 mm = r\- - (r tj ■ Vyf/w^ and 
check if ^d 2 min < a. In this case, the collision is classi- 
cally allowed. We then propagate both test particles to 
t-mim change the direction of their relative velocity v i:) - 
in a random way (thus conserving the total momentum 
and the total energy), and propagate them back to the 
original time t n . Finally, in order to take into account 
the Pauli-blocking factors in Eq. (|14j) . we calculate the 
occupation numbers f[ and /' at the new positions and 
momenta (f- = /(r^,p£) etc.) using Eq. ([5]) with Gaus- 
sians instead of delta functions, see Eq. §7). With proba- 
bility (1 — — f'j) the collision is allowed and we keep 
the new positions and momenta, otherwise the collision 
is blocked and we keep the old ones. 

We checked that the total energy is still well conserved 
when collisions are switched on: typically we find bet- 
ter than \{E)(t) - (E) (0)\/ (E)(0) ~ 1CT 5 for all times 
considered. 



D. Initialization 

Before the simulation can start, the test-particle posi- 
tions and momenta have to be initialized. Here we as- 
sume that the system is initially in equilibrium. 

A suitable equilibrium distribution is given by the dis- 
tribution function within the Thomas-Fermi or local- 
density approximation (LDA), 

/e?(r,p) = e(3 (p2/ 2m+ y T(r) _ Al) — 1 i (16) 

since it is a stationary solution of the Boltzmann equa- 
tion [22I ] . This distribution has two parameters, namely 
the inverse temperature f3 — 1/T and the chemical po- 
tential fi. The temperature T is an input parameter, 
whereas the chemical potential (i is determined by de- 
manding that the integral of Eq. (fl~6)) over r and p gives 
the right number of atoms. 

Having determined the chemical potential (i, we ran- 
domly generate the test-particle positions and momenta 
in such a way that the probability to be at position r 
and to have momentum p is proportional to / eg (r, p). In 
practice, we do this by first generating the positions ac- 
cording to the density profile obtained from Eq. @ with 




FIG. 1: Energy distribution of the atoms (divided by the 
density of states) for T/Tp = 0.2 (see text for details). The 
system consists of N — 10000 atoms with a scattering length 
a — — 0.25372/ lo (l/(kira) = —0.5). The parameters of the 
simulation are: N = 50000, w r = IMho, w p = 1.5/Zf, , and 
At = 0.02/W 



/ = f eq . Then we generate the momenta according to 

/ eq- 

E. Tests of reliability and accuracy 

In this subsection we describe two main tests we made 
to be sure that our code is reliable. Here, we assume the 
potential to be static and, as in the rest of the paper, we 
use a spherical harmonic potential 

V(v,t)=V T (v) = hm^y. (17) 

This potential defines naturally a time scale 1/oJq, a 
length scale Iho = 1 / ^/mwo , an energy scale uio, and so 
on. 

Let us consider the energy distribution of the atoms, 

In equilibrium, the distribution should be given by 
dN/dE = g(E)/(e^ E -^/ T + 1), where g(E) is the den- 
sity of states (including the degeneracy factor 2). In 
the present case of a spherical harmonic oscillator, we 
have g(E) = E 2 /ujq. In the absence of collisions, energy 
conservation automatically implies that the distribution 
stays constant, but in the presence of collisions this test 
is a non-trivial check of the Pauli blocking in the simula- 
tion. Within the test particle method, dN/dE is obtained 
by counting the test particles in energy bins. 

In Fig. [1] we show, for T/Tp = 0.2, the initial Fermi 
distribution (solid line) and the stationary distribution 
obtained in the numerical simulation after t = 30/w 
(filled circles). The agreement between the distribution 
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generated by the simulation and the initial Fermi one 
is not perfect, but satisfactory. In order to show that 
this is not a trivial result, let us see what happens if we 
switch off the Pauli blocking in the simulation of the col- 
lision term. In this case, already after a relatively short 
time ~ 3/wo , the distribution in the numerical simulation 
(empty circles) has converged to a Boltzmann distribu- 
tion with the same number of atoms and total energy 
(dashed line). So, the stability of the Fermi distribution 
in our full simulation shows clearly that Pauli blocking 
is correctly implemented. The small deviations from the 
ideal Fermi distribution are a consequence of the fact that 
with the chosen widths of the Gaussians (w r = 1.5lho and 
w p = 1.5/lho), the condition (fl"0| is not well satisfied at 
T/Tf = 0.2. When we did the same kind of compar- 
ison at higher temperatures, we found that the agree- 
ment between the simulation and the Fermi distribution 
improves: at T/Tp = 0.4, it is already perfect. 

The test described above is independent of the actual 
number of collisions. In order to check the latter, let us 
look at the collision rate 

x//i(l-/')(l-/i). (19) 

In the numerical simulation, this quantity can be ob- 
tained as Ncoii = {N/N)N co u, where N co ii denotes the 
number of collisions of test particles per unit time. 

Although in equilibrium the net effect of collisions is 
zero, the collision rate in equilibrium is a good test for 
the simulation because it can be compared with the exact 
result [Eq. (|A1|) . see Appendix For testing purposes, 
it is useful to compare also the total rate of allowed and 
blocked collisions with the exact result [Eq. (|A2[) ] . 

In Fig. [2j the collision rates (with and without block- 
ing) of the simulation are shown together with the exact 
results as functions of the temperature for two different 
values of the scattering length. In the case of relatively 
weak interaction, 1/kpa = —1 (upper panel), we see 
that the agreement between the simulation and the ex- 
act result is excellent for temperatures above ~ 0.35 Tp. 
Below that temperature, the collision rate in the simu- 
lation with Pauli blocking gradually becomes too high 
since the finite widths of the Gaussians (w r = l.blho, 
w p = 1.5/lho) do not satisfy any more the condition (TTU)) 
and act in the Pauli-blocking factors like an enhanced 
temperature. In the rest of this paper, we will therefore 
restrict ourselves to temperatures above 0.2 Tp. Near 
unitarity {1/kpa = —0.1), we consider only temperatures 
above 0.3 Tp because this is close to the supcrfluid tran- 
sition temperature at unitarity (l9j . As it can be seen in 
the lower panel of Fig. [21 the agreement between the col- 
lision rate obtained in the simulation and the exact one 
is satisfactory in the temperature range considered. The 
agreement is not as good as for 1/kpa = — 1 at high tem- 
perature because of the larger cross section which leads to 
collisions between test particles which are further apart. 
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FIG. 2: Simulated collision rates (filled circles) with and with- 
out blocking compared with the corresponding exact results, 
Eqs. (|A1|I and Eq. (jA2|) (solid lines), for a gas of 10000 atoms 
with interaction strength 1/kpa = —1.0 (top) and —0.1 (bot- 
tom). 

III. SIMULATION OF COLLECTIVE MODES 
A. Sloshing mode 

The sloshing mode is an oscillation of the center of 
mass of the system. It plays a special role because in 
a harmonic trap it is undamped and its frequency is 
equal to that of the trap, independently of the number 
of atoms, of the temperature, and of the interaction be- 
tween the atoms (Kohn mode [l?], HH). This is why it 
is often used for the experimental determination of the 
trap frequency [f|. Within the test- particle method, this 
general theorem is satisfied and it is easy to see why: 

Let us first neglect collisions. From the equations of 
motion of the individual test particles in the harmonic 
potential (fTT)) , 

r i = Pi/ m an d Pi = —mu^Ti (20) 

it is evident that the averages (r) and (p) obey analogous 
equations of motion, 

^(r) = M and l( p ) = -mu 2 {v) . (21) 
at m at 

Let us now consider the effect of a collision of two test 
particles. Of course, the trajectories of the colliding test 
particles will not obey any more the original equations of 
motion (|20[) . but the collision has absolutely no effect on 
the averages: Since the positions do not change during 
the collision, (r) remains unchanged, and since the total 
momentum of the two colliding test particles is conserved, 
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FIG. 3: Top: simulation of the sloshing mode. The mode 
was excited at t = by displacing all test particles by Iho in 
the x direction. Bottom: simulation of the breathing mode. 
The mode was excited by changing at t = all test-particle 
momenta according to p ; — > p i + cvi (c = 0.2maJo). Both 
simulations were done for a system of TV = 5000 particles at 
T = 0.4 T F and l/k F a = -0.3. 

the average (p) is not changed either. So, the equations 
of motion (|2"Tj) for the averages (r) and (p) remain valid 
in the presence of collisions. Their solution is of course 
an undamped oscillation of the center of mass (r) with 
frequency ujq. This is confirmed by the numerical result 
shown in the upper panel of Fig. [3J 

B. Breathing mode 

A couple of experiments studied the damping of the 
longitudinal and radial breathing modes [lHH 0, H[ in 
elongated traps. In a spherical trap, there is only one 
breathing mode (monopole mode), corresponding to an 
oscillation of the mean-square radius (r 2 ) around its equi- 
librium value (r 2 ) e q- In a spherical harmonic trap, this 
mode is undamped and its frequency 2ujq is independent 
of the number of collisions, like in the case of the sloshing 
mode. 

Again, this is easy to see. Consider the average ki- 
netic and potential energies, (Ekin) = (p 2 )/2m and 
(Epot) = mwg(r }/2. In equilibrium, both arc equal 
(virial theorem). Now let us assume that the system 
is compressed or expanded, such that {Eu n ) (Epot)- 
Using again the equations of motion (|20p . one obtains 

j t i(E km ) - (E pot )) = -2lu 2 q (t • p) , (22) 

j t (r-p)=2((E km )-(E pot }). (23) 

Obviously, these two equations describe an undamped 
oscillation with frequency 2luq- Let us now look if they 
stay valid in the presence of collisions. Since the collisions 
do not change the positions of the particles and conserve 
the total kinetic energy, it is clear that (Ekin) and (E pot ) 



are not affected. Now let us write the difference of (r ■ p) 
before and after a collision of two test particles i and j: 

(r • p)' - (r • p) = ir w ■ (q^ - % ) , (24) 

where qy and q^ are the relative momenta (e.g., q^ = 
(Pj — Pj)/2) before and after the collision. In the original 
collision term as written in Eq. (fT^tj) . particles have to 
be at the same position to collide, i.e., r,y = 0, such 
that (r • p) is not changed. In our simulation this is 
somewhat different, since the test particles can collide 
at a distance of up to ^Jajix. This adds a small noise 
to (r ■ p). In all practical cases, however, this noise is 
completely negligible. As an example we show in the 
lower panel of Fig. [3] the oscillation of the mean-square 
radius of the cloud function of time. As one can 
see, it is a perfectly undamped harmonic oscillation with 
frequency 2wo. 

C. Excitation of an arbitrary mode 

For the theoretical investigation of collective modes, it 
is convenient to consider a system which is in equilibrium 
until it is excited by a short pulse at t = 0. Formally, 
this is achieved by adding to the time-independent trap 
potential a perturbation term of the form 

Vx(r,t) = V 1 (v)8(t). (25) 

The reason for this choice, which is of course different 
from the experimental way of exciting a collective mode, 
is the following: Provided the perturbation V\ is small 
enough (such that the system reacts linearly to it), the 
response to a perturbation with arbitrary time depen- 
dence, Vi(r,t) = Vi(r)F(t), can easily be obtained by 
folding the result for the perturbation (f25"j) with the func- 
tion F(t). 

By integrating the Boltzmann equation over the (in- 
finitesimal) duration of the pulse, one can show that the 
effect of the perturbation (|23|) is to change the distribu- 
tion function as 

/(r,p,0 + ) = /(r,p + Vt>i(r),0-), (26) 

where + and 0~ denote the limits t — > from above and 
below, respectively. In the numerical simulation, this 
means that all test particles get a kick at t = 0, 

p i (0+)=p,(0-)-VV r i(r,(0)) (27) 

whereas their positions are not changed by the perturba- 
tion. 



D. Quadrupole mode 

From now on we will study the quadrupole mode as an 
example for a collective mode with non-trivial properties. 
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FIG. 4: Upper panel: Quadrupole response to a perturba- 
tion of the form (|28|) with c = 0.2mcjo for different tempera- 
tures. The system has TV = 10000 atoms and 1/kpa = —0.5. 
Lower panel: imaginary part of the corresponding Fourier 
transforms. 



We write the perturbation as 



y 



(28) 



corresponding to a kick at t = of p x (0 + ) = p x (0~) — 
cx(0) and p y (0 + ) = p v {0~) + cy(0). The parameter c 
determines the amplitude of the perturbation. If c is 
chosen too small, it is difficult to separate the oscillation 
of the mode from fluctuations; if it is chosen too large, 
one is not in the linear-response regime. All the following 
results were obtained with c = 0.2mujQ, corresponding to 
moderate amplitudes. By varying c within reasonable 
limits, we checked that the amplitude of the resulting 
oscillation scales linearly with c. 

After the excitation of the radial quadrupole mode, we 
can look at the time evolution of the quadrupole moment 
Q = (x 2 ) — (y 2 ) as a function of time. Results for differ- 
ent temperatures are displayed in Fig. @] Contrary to the 
sloshing and breathing modes, the quadrupole mode is 
damped and the system approaches equilibrium (Q — > 0) 
after a certain time. At high temperatures (T/Tp > 1), 
the system gets so dilute that it is in the collisionlcss 
regime {uiqt co h 3> I, t co u being the mean time between 



collisions of one atom). In this case, it takes many os- 
cillations before the system returns to equilibrium. For 
lower temperatures, the mode is damped because of the 
high collision rate (ujqt co u ~ 1), but the system is not 
yet in the hydrodynamic regime (uJt)T co u < 1) where the 
mode would become undamped again. 

For the analysis of the results, it is useful to apply a 
Fourier transform 



Q(w) 



dtQ{t)e v 



(29) 



The so-called response function is the imaginary part of 
Q{oS) and can easily be obtained from the numerical re- 
sults for Q(t) by using a fast Fourier transform (FFT) 
algorithm [26[. As an example, the Fourier transforms of 
the results discussed above are shown in the lower panel 
of Fig. 21 From the Fourier transform one can clearly 
see that the spectrum of the mode in the collisionlcss 
regime, i.e., at high temperature, has a sharp maximum 
at uj = 2wo, as it should be in an ideal Fermi gas, whereas 
at lower temperature the spectrum is broadened and the 
centroid of the spectrum is shifted to lower frequencies. 
This can be understood since at lower temperature the 
system is closer to the hydrodynamic regime, where the 
frequency should be u> = V2ujq- 

Of course, one would like to give numbers ui q and T q 
corresponding to the frequency and damping rate of the 
quadrupole mode in order to quantify these effects. The 
simplest way to obtain such numbers would be to fit the 
response function Q(t) with a damped oscillation of the 
form —Ae~ Tqt sinujqt. However, in the case of strong 
damping, this ansatz fits very badly the numerical re- 
sults for Q(t). This can be understood by looking at the 
Fourier transforms: The Fourier transform of this ansatz 
function is a Lorentzian, which has a line shape quite dif- 
ferent from that obtained in our numerical simulation for 
T/T F = 0.25 or 0.55, cf. lower panel of Fig. S Hence, 
in order to analyze our numerical results, we need some 
physically motivated ansatz for the fit. 



E. Comparison with the method of moments 

In most of the theoretical work on collective modes in 
normal-fluid Fermi gases, the Boltzmann equation was 
not solved numerically, but approximate analytical solu- 
tions were found with the help of the method of moments 
[rj li"7T - fl9| . For a detailed description of the method, see 
e.g. Ref. [H. 

Applying the method of moments to the case of a per- 
turbation of the form (|2"5)) with V\ according to Eq. (|28p , 
one obtains a theoretical prediction for the response func- 
tion ImQ(a;). A brief description of the derivation is 
given in Appendix [B] the final result reads 



ImQ(w) = — c 



8(E) 



LCT 



3m 2 (w 2 - 2w 2 ) 2 + u 2 t 2 {uj 2 ~ Alo 2 ) 2 



(30) 
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FIG. 5: Fourier transform of the numerical simulation of 
the quadrupole response (solid line) compared with Eq. (|30p 
where the parameter r is obtained by the method of moments 
(dots), or by a fit to the simulation (long dashes). The short 
dashes represent the response obtained with the extended 
method of moments including fourth-order moments. The 
system is a gas of iV = 10000 particles at T = 0.4 Tf and 
l/k F a = -0.5. 

where (E) = mwj(r 2 ) is the mean energy per atom in 
equilibrium, and r is the relaxation time as defined in 
Refs. OGIl; an d depends on the cross section (i.e., the 
interaction strength), and the equilibrium distributions, 
cf. Eq. (jBlip . One can see from Eq. ([30[) that in the col- 
lisionless and hydrodynamic limits the quadrupole mode 
has the frequencies w = 2ojq and to = v^wrj, respectively. 
The shape of the response function is completely deter- 
mined by a single parameter, r. 

By looking for the poles of Eq. ([3lJ[) . one can calcu- 
late the inverse Fourier transform which gives Q(t). The 
result has the form 

Q{t) = -Ae- r ^ smujgt + Bie^o* cosuj q t-e- rit ) , (31) 

i.e., it is a superposition of a damped oscillation with 
frequency u q and damping T qi and a non-oscillating, ex- 
ponentially decaying term. The explicit expressions for 
T\, T q , and uj q as functions of r as well as for the ampli- 
tudes A and B are given in Appendix [Cj We will refer 
to Uq and T q as the frequency and damping rate of the 
quadrupole mode. Note that in experiments determin- 
ing these quantities, the data are usually fitted with a 
function that is similar to Eq. ([31]) Q ■ 

In Fig. [5] we compare the response function obtained 
from the numerical simulation (solid line) with the result 
obtained from the method of moments, Eq. ([50]) (dotted 
line). As one can see, the height of the peak and its gen- 
eral shape are in good agreement, but the position of the 
maximum is at different frequencies. However, if we try 
to fit the numerical result with a function of the form of 
Eq. ([30]) , using r as fitting parameter, we can very well 
reproduce the numerical response function (long-dashed 
line). It is remarkable that by adjusting only one param- 
eter, t, one can simultaneously reproduce the position, 



TABLE I: Relaxation time, frequency, and damping of the 
quadrupole mode as obtained from the method of moments 
and fitting the results of the numerical simulation with a func- 
tion of the form (|30f) . corresponding to the dotted and dashed 
curves in Fig. [5] 



method 






r,/o;o 


moments 


0.451 


1.676 


0.353 


simulation 


0.587 


1.787 


0.336 



the height and the width of the peak, and also the shape 
far away from the maximum. However, surprisingly, the 
fitted value of r is larger by approximately 30% than the 
one obtained by the method of moments. As a conse- 
quence, the frequency u q and damping rate T q obtained 
from the fit of the response function deviate significantly 
from those predicted by the method of moments. These 
results are summarized in Table |TJ 

In the existing literature the method of mo- 

ments was limited to second-order moments, as described 
in appendix [Bj However, as we have seen above, this 
implies that the system is characterized by a single re- 
laxation time r, whereas in the spirit of a local-density 
approximation one would expect that in a trapped sys- 
tem the relaxation time should be position-dependent, 
t = r(r). For instance, one could imagine that the gas 
in the center of the trap is more or less hydrodynamic 
(short relaxation time), whereas far away from the trap 
center it gets very dilute and hence collisionless (long 
relaxation time). In the case of the quadrupole mode, 
this means that the Fermi-surface deformation is stronger 
at larger radii than in the trap center. It seems there- 
fore natural to include into the ansatz for the perturbed 
distribution function in addition to the standard term 
cx p x — p y describing the Fermi-surface deformation, a 
term cx r 2 (p x — p y )- More generally speaking, we should 
go beyond the standard approximation to include only 
second-order moments, and include also fourth-order (or 
perhaps even higher) moments. 

The task of extending the method of moments to the 
next higher order is in principle straight-forward but in 
practice very tedious: In the case of the quadrupole 
mode, the number of moments is increased from three 
to twelve. Some details are given in appendix [D] The re- 
sulting response function is shown in Fig. [5] as the short- 
dashed line. Surprisingly, its shape is still similar, but 
now the position of the maximum agrees rather well with 
the result of the numerical simulation (solid line). The 
agreement is even better at higher temperature (see Fig. 
[6]). This nicely confirms the correctness of our numeri- 
cal simulation and shows explicitly that the method of 
moments, if truncated at the lowest order, is insufficient. 

By doing calculations for various interaction strengths 
and temperatures, we found that the relaxation time 
from the simulation is systematically longer than that 
from the method of moments (without fourth-order mo- 
ments), Eq. ([Blip . Results for weaker and stronger in- 
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FIG. 6: Same as Fig. [5] for a temperature T = 0.7 T F . For 
clarity, the fit of the simulation is not shown. The system is 
a gas of N = 10000 particles at 1/kwa = —0.5. 



teractions {I/Ufo, = —1 and —0.1) are displayed in Fig. 
[7J We see that the general behavior of r as a function 
of temperature is the same within the simulation and 
the method of moments, but quantitatively there is a 
discrepancy of the order of 30% in the whole range of 
temperatures where our numerical simulation is very ac- 
curate (T > 0.35 Tp, cf. Fig. [5] showing the temperature 
dependence of the collision rate) . Note that at lower tem- 
peratures, the determination of the Pauli-blocking factors 
in the simulation of the collisions is not completely ac- 
curate, as discussed below Fig. [2j such that the collision 
rate below 0.35 Tp is slightly too high. Nevertheless the 
inverse relaxation time is too small. From this one can 
conclude that if we could improve the Pauli blocking in 
the simulation, the discrepancy between the simulation 
and the method of moments (without fourth-order mo- 
ments) would be even worse. The fourth order is thus 
important for the determination of the relaxation of the 
system and particularly for the frequency and the damp- 
ing of collective modes. 



F. Frequency and damping of the quadrupole mode 

As we have just seen, the numerical simulation gives 
systematically a longer relaxation time t than the 
method of moments. As the frequency uj q and damp- 
ing rate T q of the quadrupole mode are parametrized 
in terms of r (see Appendix Sec. [C|, one can ask the 
question how strongly this difference in r will affect the 
results for uj q and T q . Since we are mainly interested 
in the intermediate regime lot ~ 1 between the hydro- 
dynamic and collisionless limits, a difference of 30% in r 
can completely change the temperature dependence of ui q 
and r g . This is shown in Fig. where the crosses are the 
results obtained from the simulation, whereas the solid 
lines are the results from the method of moments. One 
can clearly see that the numerical results stay close to 
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FIG. 7: Comparison of the inverse relaxation time, 1/r, as ob- 
tained from the simulation (crosses) and from the method of 
moments Eq. (|B11[) (solid lines), as a function of temperature. 
The system consists of JV = 10000 atoms with 1/kFa = — 1 
(upper panel) and —0.1 (lower panel). 



the collisionless limit to much lower temperatures than 
the results obtained by the method of moments. 

To estimate the resulting precision of our numerical 
result on ui q and T g , we show in Fig. [8] the error bands 
(short-dashed lines) which we obtain if we assume that 
our simulation may give a r which is wrong by at most 
15%. This error includes numerical uncertainties which 
can be estimated from the scattering of the points in Fig. 
[7] and the systematic deviation of the collision rate shown 
in Fig. H 

If we include the fourth order moments, a global relax- 
ation time t does not exist anymore but we could define 
an effective one by fitting the response function as we did 
for the simulation. This effective relaxation time agrees 
very well with the one of the simulation so that both re- 
sults give very similar frequency and damping. However, 
from a theoretical point of view, the definition of these 
quantities should come from the zeroes of the determi- 
nant of the matrix Ay defined in Appendix [Dj Such a 
discussion is postponed to a forthcoming publication [29j ] . 
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FIG. 8: Frequency (top) and damping rate (bottom) of the 
quadrupole mode as a function of the temperature as obtained 
from the numerical simulation (crosses) and from the method 
of moments (solid line). The short-dashed lines indicate the 
error band of the results from the simulation if we admit that 
the relaxation time r of the simulation may be wrong by 15%. 
The system consists of N = 10000 particles close to unitarity 
{l/k F a = -0.1). 



IV. CONCLUSIONS 

In this paper, we presented a test-particle method for 
solving numerically the Boltzmann equation for trapped 
Fermi gases. While such methods have been popular in 
other fields of physics for many years, there have been 
only a few applications to ultracold atomic gases I14J- 

HE 

l2ll |. Our method is similar to that of Refs. [Til Il6l| 
with some differences in the treatment of the collision 
term. In order to compute the occupation numbers in the 
Pauli-blocking factors in the collision term, we represent 
each test particle by a Gaussian in r and p space. The 
minimum value of the width of the Gaussian is dictated 
by the statistical fluctuations due to the finite number of 
test particles, limiting the applicability of the method to 
temperatures above ~ 0.2 Tp. 

As a first application of the method we discussed some 



collective modes. For simplicity, we considered only toy 
systems consisting of ~ 10 4 atoms in a spherical har- 
monic trap and neglected the mean-field potential and 
medium modifications of the cross section. As expected, 
the sloshing and monopole modes are undamped and in- 
dependent of the collisions. In contrast, the quadrupole 
mode is very sensitive to collisions. In the hydrodynamic 
limit, its frequency should approach y/2u>o, while it is 
2wo in the collisionlcss limit. In our simulations, we 
never reach the hydrodynamic regime, but the collision- 
less regime can be realized at high temperature due to 
the diluteness of the gas. 

Surprisingly, the frequency and damping rate of the 
quadrupole mode obtained within the numerical simula- 
tion are quite different from those obtained within the 
widely used method of moments including moments up 
to second order in r and p. The method of moments 
predicts a relaxation time r which is significantly shorter 
than the one obtained within the simulation. The rea- 
son is that the r dependence of the relaxation time is 
neglected if only the p\ — p y moment is taken into ac- 
count for the description of the Fermi-surface deforma- 
tion. We have shown that if the method of moments is 
extended to moments up to fourth order in r and p, e.g., 
the r 2 (p 2 — Py) moment, the agreement with the simula- 
tion becomes very good. 

The focus of the present paper was mainly to explain 
the test-particle method and to show its usefulness. For 
instance, the deficiency of the method of moments up to 
second order would not have been detected without the 
comparison with the numerical result. In future studies, 
we plan to apply the method to more realistic cases. In 
particular, in order to reach the typical numbers of atoms 
in the experiments, we will have to increase A by a factor 
of ~ 10 — 100. However, this should not pose a big prob- 
lem: According to Eq. (jTU)) , if we increase TV but keep the 
ratio T/Tp fixed, the widths w r and w p may be chosen 
larger (oc A 1 / 6 ). This means that Eq. §§§ stays satis- 
fied with the same number of test particles, A, i.e., with 
a reduced ratio A/A. The computation time will only 
grow because of the increased collision rate (due to the 
larger test-particle cross section a = aN/2N). Another 
point is the trap geometry. The traps in the experiments 
are usually not spherical, but elongated. Concerning the 
propagation of the test particles, this does not cause any 
difficulty, but in the calculation of the occupation num- 
bers, it will probably be necessary to replace the width 
w r of the Gaussian in r space by different widths w x , w y , 
and w z in the three space directions. Another impor- 
tant advantage of the numerical method is that an an- 
harmonicity of the trap potential, which is always present 
in real experiments, can easily be included. 

Finally, the mean field II 9i and medium modifications 
of the cross section [|| Il9l | should be included. The 
mean field, which originally depends on the chemical 
potential /i and the temperature T, can be expressed 
as a function of the local density and energy density, 
which are both obtainable in the simulation. However, 
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as shown previously [lj|, the mean field is not just pro- 
portional to the density: this leads to a huge numerical 
effort which is beyond the scope of this paper. The in- 
medium cross section is also difficult to be included be- 
cause it depends on too many variables to be tabulated: 
(7 = a{k = |p + Pil/2,3 = |p - p x |;/U,T). One possi- 
ble solution of this problem is to replace the full k and 
q dependence of the in-medium cross section by a simple 
paramctrization which results in the same local relax- 
ation time r(/i,T). Work in this direction is already in 
progress. In Refs. [1, [HI, within the method of moments 
up to second order, the use of the in-medium cross section 
spoiled the agreement with experimental data because 
the resulting relaxation times were too short. Since the 
present work shows that the numerical simulation gives a 
longer relaxation time than the method of moments, we 
hope that this problem can be solved. 

Further important extensions of the present work are 
the generalization to polarized Fermi gases and to super- 
fluid systems. These questions, however, require more 
fundamental theoretical studies before they can be tack- 
led numerically. 

Appendix A: Collision rate at equilibrium 

Replacing the distribution functions in Eq. (fTT)| by 
equilibrium distribution functions f eq , one obtains after 
some algebra the equilibrium collision rate 



coll,eq 



1 

i^ 1 



d 3 r 



dkk 2 I dqq 2 —a(q) 
o Jo m 

/ tanh~ (tanh =j tanh y) \ 



YsinhX 



(Al) 



where k = p + p 1 , q = (p p x )/2, X = /3(k 2 /8m + 
q 2 /2m+ Vr — (J,) and Y = fikq/2m. The total rate of al- 
lowed and blocked collisions is, in turn, given by Eq. (fH))) 
but without the factor (1 — /')(1 — /{) in the integrand, 
leading to 



7y"C+ blocked) 
coll,eq 



1 

4^ 



d 3 r [ dkk 2 f dqq 2 --a(q) 
Jq Jo m 

tanh _1 (tanh =y tanh -y) 
Ye x sinh X 



(A2) 



If the trap potential is spherically symmetric or har- 
monic, the spatial integrals can be reduced to one- 
dimensional ones. The remaining three-dimensional in- 
tegrals are evaluated numerically with a Monte-Carlo al- 
gorithm. 



one in the form 

/ - feq = feq(l - feq)$(r,p,t) 



(Bl) 



Inserting this expression into the Boltzmann equation 
and keeping only terms linear in the perturbation, one 
obtains (see Eq. (36) of Ref. for the case without 
mean field but with an external perturbation): 



f eq (l - f eq ) $ + ^ • V r $ - V r V T • V p $ 

+ J gZ.V r 7i) =-/[*]• (B2) 
m / 

Here, /[$] is the linearized collision term as defined in 
Eq. (37) of Ref. [l9[ (up to a factor (2tt) 3 since here we 
are using a different normalization of /): 



d 3 pi 
(2tt) 3 



e^^|v - Vi|/ eg / eg i 



The perturbation Vi is given by Eqs. ([25J and (|28|). The 
usual approximation consists in making the ansatz 



x (1-/')(W' !)(* + *i -$'-*!). (B3) 



*(r,p,t) = 5^c i (t)0i(r I p). 



(B4) 



i=l 



with time-dependent coefficients Cj and <f>i — x 2 — 
y 2 , <j> 2 = xp x ~ VPy, and (j> 3 = p 2 x - p\, i.e., only 
quadratic moments are considered. Evaluating the mo- 
ments / d 3 rd 3 p(f>i(r, p) x Eq. (|B2|) . one obtains a system 
of equations for the Fourier transformed coefficients: 



(B5) 



with 



d 3 rd 3 p 
(2tt) 3 ' 

iui(j)j - 



feq (1 feq ) 



31 2m 



9 1 



(B6) 



and 



d 3 rd 3 p 
(2tt) 3 ' 



ifeq(l ~ feq)~ ■ W X , (B7) 



where {.,.} are the Poisson brackets. Using the virial 
theorem, we obtain explicitly: 



Appendix B: Quadrupole response within the 
method of moments 



In the case of a weak perturbation, we can write the de- 
viation of the distribution function from the equilibrium 



-icoci — muj C2 = 



2ci — imwc2 — 2m 2 u>QC3 



-13c, 



C2 



iui 



(B8) 
(B9) 

(BIO) 
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where the relaxation time r is defined by @, G3 



1 



3/3 



r m 2 A(£ fcm ) J (2tt) 
Solving this system of equations, we find 
f3uj 2 c(l — iijjt) 



ci(w) = 



oj 2 — 2cJq — iu!t(uj 2 — 4Wq) 



(B12) 



and similar expressions for c 2 and C3. However, only the 
coefficient c\ contributes to Q: With Eq. (IB1|) . and using 
again the virial theorem, we obtain Q(t) = (x 2 — y 2 ) = 
4:T(r 2 )ci(t)/3muj 2 , or explicitly: 



Q(u) 



4(r 2 )c 



1 — IU1T 



3m ui 2 — 2u> 2 — iwr(w 2 — 4,lo 2 ) 
Taking the imaginary part, we obtain Eq. ([30 



(B13) 



Appendix C: Time dependence of the quadrupole 
response within the method of moments 

In order to compute the Fourier transform of Eq. 
(|B13|) . let us start by factorizing the denominator: 

Ld 2 — 2co 2 — iu!T(uj 2 — 4cJq) 

= —ir(ui — lui)(uj — lu 2 )(uj — ^3) • (CI) 

The expressions for the roots Wj can be given in closed 
form. Defining f = ujqt and 



6 = ( 1 + 9f 2 + 3f - 39-r 2 + 192-r 



1/3 



1 / 1 - 12f 2 



(C2) 
(C3) 



we can write the roots LOi as 

u>i = -iTi , u 2 = Lo q - iV q , u> 3 = -u q - iT q , (C4) 
with 

r i = T^I + u + , r <? = - -y > w ? = — ■ ( C5 ) 



Now it is straight-forward to evaluate the inverse Fourier 
transform of the response (|30[) using the residue theorem. 
The result is given by Eq. ([31]) with 



.4 



4 C (r 2 ) ^r + crr-iYKi-ryr) 



3mw 9 T 
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u 2 + (r, - ro 2 

4c(r 2 ) 1 - r x r 



3mr ^ 2 + (r g -r 1 ) 5 



(C6) 
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Appendix D: Extension of the method of moments 
to fourth-order moments 



Taking fourth order moments into account, we extend 
the previous ansatz Eq. (|B4[) as follows: 



ciO 2 - V 2 ) + c 2 (xp x -yp y ) + c 3 {p 2 x 



Pi) 



+ c A r 2 {x 2 - y 2 ) + c 5 p 2 (x 2 - y 2 ) + c 6 r ■ p(x 2 
+ c 7 r 2 (xp x - yp v ) + c 8 p 2 (xp x - yp y ) 
+ c 9 r • p(xp x - yp y ) + c 10 r 2 (p 2 x - p 2 y ) 



+ c llP 2 (pl - pi) + c 12 r ■ p(p 2 x 



Pi) 



y 2 ) 



(Dl) 



which can be written as $ = J2i=i Cifii with, for example, 
4>i = (x 2 — y 2 ). Following the same steps as explained in 
Appendix[B] we obtain now a system of twelve equations. 
The matrix A^j can be computed explicitly. Contrary to 
the second order calculations, the virial theorem can no 
longer be used to reduce the number of unknown quan- 
tities so that the system now depends on (r 2 ), (r 4 ) and 
(r 6 ). In the matrix elements of the collision term, more 
parameters appear, generalizing the single parameter r 
of the second order method. 

After solving the system of equations numerically, we 
can express the quadrupole moment in terms of the co- 
efficients Ci as : 



Q(w) 



4T 
~3~ 



7c 4 



3rac5 + meg 



(D2) 

Further details and explicit formula for the matrix can 
be found in (H. 
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